function cfdTimeAdvanceLUSGSUp
%--------------------------------------------------------------------------
%  Written by LiuHaHa @ NWPU, 2022 01 
%  Contact me at: liuzhikan@mail.nwpu.edu.cn
%==========================================================================
% Routine Description:
%   This function calculate Lower LU-SGS element
%--------------------------------------------------------------------------
global cfdFlow;

theNumberOfElements = cfdFlow.mesh.numberOfElements;
theneighEleNum = cfdFlow.mesh.neighEleNum;
theLUSGSEleLabel = cfdFlow.mesh.LUSGSEleLabel;
theneighEleIndex = cfdFlow.mesh.neighEleIndex;
theneighFaceIndex = cfdFlow.mesh.neighFaceIndex;
cfdBoundaries = cfdFlow.mesh.cfdBoundaryPatchesArray;
numberOfBoundaries = cfdFlow.mesh.numberOfBoundaryPatches;
thenumberOfInteriorFaces = cfdFlow.mesh.numberOfInteriorFaces;
gama = cfdFlow.gama;
owners = cfdFlow.mesh.owners;
faceSf = cfdFlow.mesh.faceSf;
eign_f = cfdFlow.eign_f;
eignV = cfdFlow.eignV;
phyTemp = cfdFlow.Phy;
dltPhy = cfdFlow.dltPhy;

for iElement=theNumberOfElements:-1:1
%     tempFlux = zeros(slvpara.ND+2,1);
    dltPhy(iElement,:) = dltPhy(iElement,:)*eignV(iElement,1);
    for i=1:theneighEleNum(iElement,1)
        if(theLUSGSEleLabel(iElement,i) == 1)
            iNeigFace = theneighFaceIndex(iElement,i);
            iNeigEle  = theneighEleIndex(iElement,i);

            phyPre_f = phyTemp(iNeigEle,:);
            phyNxt_f = phyPre_f + dltPhy(iNeigEle,:);

            PressurePre = (gama-1)*(phyPre_f(end) - 0.5*dot(phyPre_f(2:end-1),phyPre_f(2:end-1))/phyPre_f(1));
            PressureNxt = (gama-1)*(phyNxt_f(end) - 0.5*dot(phyNxt_f(2:end-1),phyNxt_f(2:end-1))/phyNxt_f(1));

            Vn_pre = dot(phyPre_f (2:end-1),faceSf(iNeigFace,:))/ phyPre_f (1);
            Vn_Nxt = dot(phyNxt_f (2:end-1),faceSf(iNeigFace,:))/ phyNxt_f (1);
            if(owners(iNeigFace) == iNeigEle)
                Vn_pre = -Vn_pre;
                Vn_Nxt = -Vn_Nxt;
            end

            dltFlux = -0.5*(phyNxt_f*Vn_Nxt - phyPre_f*Vn_pre);
            if(owners(iNeigFace) == iNeigEle)
                dltFlux(2:end-1) = dltFlux(2:end-1) - 0.5*(PressureNxt - PressurePre)*(-faceSf(iNeigFace,:));
            else
                dltFlux(2:end-1) = dltFlux(2:end-1) - 0.5*(PressureNxt - PressurePre)*(faceSf(iNeigFace,:));
            end
            dltFlux(end) = dltFlux(end) - 0.5*(PressureNxt*Vn_Nxt - PressurePre*Vn_pre);

            dltPhy(iElement,:) = dltPhy(iElement,:) + dltFlux(:)';
            dltPhy(iElement,:) = dltPhy(iElement,:) + 0.5*eign_f(iNeigFace)*dltPhy(iNeigEle,:);
        end

    end
%     dltPhy(iElement,:) = dltPhy(iElement,:) + tempFlux(:,1)';
    if( eignV(iElement,1) > 1.0e-15)
        dltPhy(iElement,:) = dltPhy(iElement,:)/eignV(iElement,1);
    else
        dltPhy(iElement,:) = 0.0;
    end
end
for iBPatch=1:numberOfBoundaries
    theNumberFacesOfBoundaryPatch = cfdBoundaries{iBPatch}.numberOfBFaces ;  
    startFace = cfdBoundaries{iBPatch}.startFaceIndex;
    for iBface=startFace:startFace+theNumberFacesOfBoundaryPatch-1
        iBelement = iBface - thenumberOfInteriorFaces + theNumberOfElements;
        dltPhy(iBelement,:) = 0.0;
    end
end

cfdFlow.dltPhy = dltPhy;













